
GRAVITATIONAL EXPERIMENTS 
WITH A COLLISIONLESS 
TWO-DIMENSIONAL COMPUTER MODEL 


by Frank Hohl and Stephen K. Park 

Langley Research Center 
Langley Station, Hampton, Va. 



NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


WASHINGTON, D. C. ♦ JULY 1968 


1 


TECH LIBRARY KAFB, NM 




TECH LIBRARY KAFB, NM 



□131320 


GRAVITATIONAL EXPERIMENTS WITH A COLLISIONLESS 

TWO-DIMENSIONAL COMPUTER MODEL 

By Frank Hohl and Stephen K. Park 

Langley Research Center 
Langley Station, Hampton, Va. 


NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 


For sale by the Clearinghouse for Federal Scientific and Technical Information 
Springfield, Virginia 22151 — CFSTI price $3.00 



GRAVITATIONAL EXPERIMENTS WITH A COLLISIONLESS 
TWO-DIMENSIONAL COMPUTER MODEL 

By Frank Hohl and Stephen K. Park 
Langley Research Center 

SUMMARY 

A two-dimensional model is used to perform computer experiments for a collision- 
less self -gravitating system. The gravitational field is obtained by solving the Poisson 
equation and the system is advanced stepwise in time. Computer simulations have been 
performed for systems with an initially uniform distribution over a circular region in 
x,y space, zero thermal velocity, and various values of initial solid-body rotation. Up 
to 4000 stars were used in the calculations. 

INTRODUCTION 

Actual stellar systems, such as galaxies, contain about loH stars with a sufficiently 
large average separation so that binary encounters between stars can be neglected. Stel- 
lar systems can therefore be described by the collisionless Boltzmann equation. In order 
to simulate the evolution of such stellar systems on a computer, the motion of at least 
several thousand masses or stars should be followed. Hohl and Feix (refs. 1 and 2) and 
Lecar and Cohen (ref. 3) used one -dimensional sheet models to study the evolution of self- 
gravitating systems. A similar model where the motion of a large number of concentric 
spherical mass shells is followed has been used by Henon (refs. 4 and 5). Recently, 
Hockney* (ref. 6) and Hohl (refs. 7 and 8) introduced a two-dimensional model where the 
stars are represented by infinitely long mass rods. In the present report the results of 
some simple experiments with the two-dimensional model are given. The calculations 
show that filamentary and other irregular structures similar to those of some actual 
stellar systems can be obtained by purely gravitational effects. 

SYMBOLS 

ajj variable defined by equation (3) 

D impact parameter 

*The results presented in the present report and in reference 6 are similar and they 
were obtained independently by the present authors and by Hockney. The methods used in 
the numerical calculation are different. The present method is more general and much 
simpler to use than the Fourier method used by Hockney. 
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gravitational field, -Vcp 
summation index 
mass per unit length 
star density 

number of stars in system 
potential energy 
magnitude of radius vector 
radius vectors defining position 
real part 
time 

kinetic energy 
total energy 
velocity 

transverse velocity 
positions along X- and Y-axes 
complex position, x + iy 
overrelaxation parameter 
variable defined by equation (15) 
mass density 
relaxation time 
equilibrium rotation period 
gravitational potential 





frequency of rotation 


u)g equilibrium frequency of rotation 

Subscripts: 

j,k,n,m summation indexes 

max maximum 

min minimum 

r radial 

x,y x- and y-component 

9 azimuthal 

Arrows over symbols denote vectors. < > denotes expectation values. 

DESCRIPTION OF MODEL 

The model consists of a large number of rods of equal mass per unit length and the 
rods are of infinite extent in the z -direction. These rods move in the x,y plane under 
the action of their mutual gravitational attraction. The system of mass rods is advanced 
in time in the following manner. First, the distribution of mass p(x,y) is used to obtain 
the gravitational potential cp(x,y) by numerically solving the Poisson equation. Second, 
the gravitational field at the position of the particles is computed from the potential 
cp(x, y). Third, Newton’s laws are used to advance the motion of all the mass rods for a 
small time step 6t. These three steps represent one cycle and they are repeated until 
the desired evolution of the system is achieved. 

The crucial point in the computations is the solution of the Poisson equation. It is 
desirable that the time required for this process be only a small fraction of the cycle 
time. If the system is advanced for a small time step 5t, the mass distribution p(x,y) 
will not change very much. The change in the gravitational potential will then also be 
very small. Thus, the solution of the finite difference form of the Poisson equation by an 
iteration method which uses the potential from the previous cycle as an initial guess will 
converge very rapidly. It was found that 5 to 7 iterations per cycle gave satisfactory 
results. 

To solve the Poisson equation, the boundary conditions around the rectangular mesh 
used in solving the finite difference form of the Poisson equation are required. The 
potential at an arbitrary boundary point z = x + iy can be written as 
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(la) 


N 

cp{K,y) = 2Gm ^ log e |z - z n | 
n=l 


where z n = x n + iy n is the position of a rod of mass m, G is the gravitational con- 
stant, and N is the number of particles in the system. 

Equation (la) can be written as 


<p(x,y) = 2GmN log e |z| + 2Gm Re 


I loge (l - ) 

n=l 


Since 



the potential of the boundary points can be approximated by 


(lb) 


where 


<p(x,y) = 2GmN log e |z| - 2Gm Re 



a k 

kz^ 


a k = 



( 2 ) 


(3) 


and the series expansion for log e ^1 - ~^j has been truncated after 15 terms. The 

accuracy of equation (2) has been checked by direct summation of the log e |z - z n | 
potential and was found to agree to better than 0.1 percent. Equation (2) for the poten- 
tial on the boundaries has also been used by Hockney (ref. 6). Equation (la) can also 
be used to obtain the potential of an arbitrary point. However, the time required for 
such a method is too large. Therefore, the potential is obtained by solving the Poisson 
equation. 

The Poisson equation 



+ = 4?rGp(x,y) 

8y^ 


( 4 ) 


is solved by using the standard five-point difference equation (ref. 9) 
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^n+ljin + ^n,m+l + ^n-l,m + ~ ^GP^m 


(5) 


where Ax = Ay = 1. The mass per unit length of all the rods in the cell n,m is repre- 
sented by p n m . This set of simultaneous equations is solved by an iteration method on 
the Control Data 6600 computer system in the form 

^n,m = ^n,m + ^^°n-l,m + ^n+ljm + < Pn,m-l + Vn^m+l ~ 4<Pn,m “ 47rGp n>m ^ (®) 

To save computer storage and increase the convergence rate, the new values of cp r+ 1 
which have already been determined are used in the right-hand side of equation (6). The 
superscript r refers to the rth iteration and the parameter y is adjusted to give the 
maximum rate of convergence. 

The differential equations of motion for the stars are 

d 2 Rj - 

-zd-HW) m 


dRj 

dt 


V, 


(8) 


where R is the position of a star and K is the gravitational field. The evolution in 
time of the particle trajectories is given by the finite -difference approximation 


2 

Rj(t + 6t) = Rj(t) + ^j(t)6t + K(Rj,t) 

(9) 

2 — 

(10) 


where 


was approximated by the backward difference 


dK _ 
dt 



(ID 


Equations which were found to be nearly as accurate as equations (9) and (10) but 
which require fewer computer operations are (from ref. 10) 
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Vj(t + 5t) = Vj(t) + 6tK(Rj,t) 

(12) 

Rj(t + 6t) = Rj(t) + 6tVj(t + 5t) 

(13) 


The term ^ in equation (10) was introduced to remove a computational instability. 
Similarly, the use of the new velocity in equation (13) removes an unconditional computa- 
tional instability. 

Equations (12) and (13) can be expressed in a simpler form by introducing the 
quantity 


?j(t) = Rj(t) - Rj(t - 5t) 


(14) 


The equations can then be written in the simplified form 

|j(t + St) = £.(t) + K(R^ St 2 


(15) 


and 


Rj(t + St) = Rj(t) + fj(t + St) (16) 

Equations (15) and (16) were also used by Hockney (ref. 6) and they require the least 
number of computer operations. However, if it is desired to display the evolution of the 
system in velocity space or to compute the kinetic energy of the system, the equation 


V(t) = [l(t) +l(t + 6t)j 


(17) 


must be evaluated. It may then be more desirable to use equations (12) and (13). All 
three sets of equations were used and all gave similar results. 


The kinetic energy of the system is 



j=l 


(18) 


and the potential energy is 
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N 

p 4 m I^i) (19) 

i=i 

where the summation goes to N, the number of stars in the system. The total energy of 
the system U remains constant and is 


U = T + P (20) 

For the calculations presented here, the Poisson equation was solved on a 51 by 
51 mesh. During the evolution of the system, the potential at several mesh points was 
computed by summing directly the contribution to the potential from each mass rod. 

These potentials agreed with those obtained from a solution of the Poisson equation to 
within 0.1 percent. 

Another check on the method of solution was made by comparing the evolution of 
two systems with identical initial conditions where the mesh size used in the solution of 
the Poisson equation was 51 by 51 for one system and 101 by 101 for the other. The 
resulting evolution of the two systems was nearly identical. The cycle time for a 
2000-particle system with a 51 by 51 mesh and with 7 iterations per cycle is 1 second on 
the Control Data 6600 computer system at the Langley Research Center. This time 
includes such operations as the checking of the potential, the calculation of energy and 
angular momentum, and writing the positions and velocities of all stars on tape. 

The relaxation time (ref. 11) for the two-dimensional model is of interest to deter- 
mine the number of rotations for which the system can be treated as collisionless. Con- 
sider an encounter between a rod of mass m per unit length and velocity V with a 
stationary rod of mass M per unit length. If D is the distance of closest approach 
(impact parameter), the transverse force acting on the moving rod during a time 2D/V 
can be approximated as 2GmM/D. (See ref. 12.) The moving star has therefore acquired 
a transverse velocity 


6V ± «^M (21) 

For a star interacting with a system containing many stars, the effect of many individual 
encounters must be summed. The number of encounters in a time t with impact param- 
eter between D and D + dD is tVn dD where n is the density of stars in the sys- 
tem. The expectation value of Vj_^ is then given by 
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= 16tnG 2 M 2 D W y Dmin 


= 16tnG 2 M 2 


( 22 ) 


since, in general, 
defined as the time 
Thus, for <Vj_J> 
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max 


» D 


min 


required for 
= V 2 , t = r r 


where 

<Vj.2> 

and 


D m ax = Ax. The relaxation time t c is 
to be of the same order of magnitude as V 2 . 


T c * 


8G 2 m 2 nD 2 iax 


(23) 


where M = m and D max is taken to be equal to the dimension of the system, which 
for systems near equilibrium is the Jeans or Debye length (ref. 2). 

The gravitational field inside a system with a uniform circular distribution in coor- 
dinate space is 


K = 2irGpR 

Balancing the gravitational force toward the center of the system by the centrifugal force 
co 2 R requires a frequency 


u> = ^2nGp 

The rotational period for such a system then becomes 


r - 2 E 

S /2irGp 

For the systems investigated, the ratio of the relaxation time r c to the rotation 
period t s is 
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Therefore, only for a time less than about 4 rotations can the system be considered as 
collisionless. To keep the effects of collisions negligible for the first 10 rotations 
requires a system containing about 10 000 stars. 

For the 2000-particle systems investigated, the total energy was conserved to 
within 2 percent. The error in the energy for systems containing 4000 particles was 
less than 0.5 percent. 


RESULTS AND DISCUSSION 

The computer experiments were performed for systems which have an initially 
uniform circular distribution in x,y space and zero thermal velocity. The evolution of 
such systems is then studied for various values of initial solid-body rotation. The initial 
conditions are obtained by using a random-number generator which gives a nearly uni- 
form distribution over a circular region of the x,y plane. The normalizations 47 tG = 1 
and m = 1 were used for all the calculations. 

First, the results for the case where the system is in equilibrium are presented. 
For this case the initial frequency of rotation u> equals u)g, the frequency required 
such that the centrifugal force balances the gravitational attraction towards the center of 
the system. Thus u> = u>g, where 


u>g^ = 27rGp 

Figure 1 shows the evolution of a 2000-particle system in x,y coordinate space. The 
time t has been normalized to the period of rotation Tg = The time step used in 
the calculations is 6t = Lamb (ref. 13) has given the result that an infinitely long 

circular cylinder of uniform density may rotate in relative equilibrium about its longi- 
tudinal axis with an angular velocity given by u>% = uG p. One would therefore expect the 
initial distribution shown in figure 1 to be unstable. Thus, at t = 0. 5Tg, there appears 
a fourth harmonic perturbation around the periphery of the system. Later the perturba- 
tion goes to an egg shape and it finally disappears and leaves the system in a steady 
state. At t = 0.75rg and t = l.OOTg, there is an indication of the appearance of spiral 
arms; however, it appears that collisional effects suppress this tendency. The evolution 
of the corresponding velocity distribution is shown in figure 2. There is almost no 
change in the velocity distribution as the system evolves in time. Figure 3 shows the 
velocity distribution of the system after the initial solid-body rotation has been sub- 
tracted. The radial velocity V r is plotted against Vg - rcUg where Vg is the azi- 
muthal velocity and r is the radius from the center of the system to a star. Initially, 
the thermal (or random) velocity of all mass rods is zero. Because of the small 
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perturbation caused by the random initial position of the stars, a small thermal velocity 
builds up and the system expands in V r - - rwg) space. After a time t = 2.0T g , 

there is little further change in the velocity distribution. The increase in the thermal 
velocity causes the system to stabilize. 

These calculations with o> = Wg were repeated for a 4000-particle system. The 
results are shown in figure 4. It can be seen that the deviation from a circular distribu- 
tion is even smaller now than for the 2000-particle system. 

The results for the case of zero initial rotation, u> = 0, are presented next. Fig- 
ure 5 shows the evolution of the system in coordinate space. After an initial implosion, 
the system expands and presents some highly filamentary structure. At t = 0.53Tg, sev- 
eral clusters of stars have condensed. However, at a later time only a large central 
cluster remains and the system takes on an appearance reminiscent of an elliptical gal- 
axy. The corresponding evolution in velocity space is shown in figure 6. Initially, all 
stars have zero velocity. As the system collapses under the gravitational forces, the 
velocity of the stars increases. The system then expands rapidly in velocity space and 
shows filamentary structure similar to that appearing in coordinate space. A motion 
picture was prepared which shows the evolution of the system. One of the striking fea- 
tures that can be observed in the motion picture is the appearance of streamers in veloc- 
ity space. Chains of stars move out of the main body of the system and then curve back 
again. After several pulsations, the pressure due to the increased temperature builds 
up and the system approaches an equilibrium state. 

For a system with an initial rotation equal to half that required to balance gravita- 
tion, the system again contracts initially and then expands. The results are shown in 
figure 7(a). At time t = 0.53Tg, the system has developed an irregular structure which, 
however, disappears again and after only a few rotations the system approaches an equi- 
librium configuration similar to that of an elliptical galaxy. Figure 7(b) shows the corre- 
sponding evolution in velocity space. 

In figure 8(a) the results are shown for a system with an initial rotation equal to 
1.3o)g. The general behavior is very similar to that of the previous case with u> = 0.5a>g. 
At the time t = 0.75Tg the system shows a structure similar to that of the Crab Nebula. 
The evolution in velocity space for this system is shown in figure 8(b). 

The time development of the total kinetic energy of the 2000-particle systems inves- 
tigated is shown in figure 9. Only small oscillations occur for the balanced case. For 
the other cases the initially large oscillations in the kinetic energy are quickly damped 
as the systems approach an equilibrium state. The rapid damping of the oscillations in 
the kinetic energy is another indication that collisional effects are important even after 
only 3 to 4 rotations. 
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CONCLUDING REMARKS 


The two-dimensional model used to study the evolution of self-gravitating systems 
indicates that a variety of filamentary and other structures can be produced by purely 
gravitational effects. Some of the systems investigated showed a structure similar to 
that of the Crab Nebula. 

The evolution of the balanced cylinder showed that there was a tendency for spiral 
structure to appear and also a tendency for the cylinder to become unstable. It appears 
that both of these effects are quickly suppressed by collisional effects. Thus, systems 
with more stars should be used to increase the ratio of the relaxation time to the rota- 
tional period. For example, to keep the system collisionless for 10 rotations requires 
a system containing 10 000 mass rods. 

The method used for the calculations in this report is general and requires no mod- 
ification either to increase the number of stars in the system or to increase the number 
of mesh points used in the numerical solution of the Poisson equation. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., April 16, 1968, 

129-02-01-01-23. 
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Figure 3.- Buildup of the thermal velocity for a cylindrical stellar system with ui = «g. 
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Figure 5.- Evolution in coordinate space of a cylindrical stellar system with cd = 0. 
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Figure 8.- Evolution of a 2000-star system with u) = 1.3o>g. 
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